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We present an approach to the implementation ol the virtual crystal approximation (VCA) for the 
study of properties of solid solutions in the context of density-functional methods. Our approach 
can easily be applied to any type of pseudopotential, and also has the advantage that it can be used 
to obtain estimates of the atomic forces that would arise if the real atoms were present, thus giving 
insight into the expected displacements in the real alloy. We have applied this VCA technique within 
the Vanderbilt ultrasoft-pseudopotential scheme to predict dielectric and piezoelectric properties of 
the Pb(Zro.5Tio.s)03 solid solution in its paraelectric and ferroelectric phases, respectively. Compar- 
ison with calculations performed on ordered alloy supercells and with data on parents compounds 
demonstrates the adequacy of using the VCA for perovskite solid solutions. In particular, the VCA 
approach reproduces the anomalous Born effective charges and the large value of the piezoelectric 
coefficients. 
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I. INTRODUCTION 

The application of first-principles electronic band- 
structure methods to the study of disordered alloys and 
solid solutions requires some approximation for the treat- 
ment of the alloy disorder. A "direct" approach is to 
make use of the supercell approximation, i.e., to study 
one or more disordered configurations in a supercell with 
artificially imposed periodic boundary conditions. Such 
calculations generally require the use of very large super- 
cells in order to mimic the distribution of local chemical 
environments, and tend to be computationally very de- 
manding. A much simpler and computationally less ex- 
pensive approach is to employ the virtual crystal approxi- 
mation (VCA) fjj , in which one studies a crystal with the 
primitive periodicity, but composed of fictitious "virtual" 
atoms that interpolate between the behavior of the atoms 
in the parent compounds. This technique has seen wide 
use in band-structure calculations P-jO]]. Another pos- 
sible approach would be to make use of the coherent po- 
tential approximation (CPA) |l2|, but unfortunately the 
CPA is generally not well suited for use in first-principles 
total-energy methods. A different way to go beyond the 
VCA is to carry out a systematic perturbation expansion 
in the difference between the true and VCA potentials, 
an approach that is sometimes referred to as "compu- 
tational alchemy" 0-0]. However, this method is much 
more complicated than the usual VCA, requiring the use 
of density-functional linear-response techniques. 

Clearly the VCA has the advantages of simplicity and 
computational efficiency, if two possible concerns can be 
addressed. First and foremost is the question of the ac- 
curacy of the VCA approximation. Previous work has 



demonstrated good accuracy for the VCA in some semi- 
conductor and ferromagnetic materials ^-0], but it was 
found to be inadequate for an accurate treatment of the 
electronic structure of some unusual semiconductor sys- 
tems J8| |l0| . Until the recent pioneering work of Ramer 
and Rappe [ fLl| , nothing was known about the ability 
of the VCA to describe the properties of an important 
class of materials, the ferroelectric perovskite solid so- 
lutions. Their work strongly suggests that these alloys 
are good candidates for modeling with the VCA, since it 
reproduces the strain-induced transitions of ordered su- 
percells of Pb(Zr .5Ti .5)O3. However, it is not known 
whether the VCA is good enough to predict the anoma- 
lous dielectric and piezoelectric properties of perovskite 
solid solutions. 

A second concern is more technical. By its nature, 
the VCA is closely tied to the pseudopotential approx- 
imation. Indeed, unless pseudopotentials are used, it is 
hopeless to apply the VCA to the usual case of isoelec- 
tronic substitution (i.e., atoms belonging to the same col- 
umn but different rows of the Periodic Table). However, 
as pseudopotential methods have matured, it has become 
less obvious what is the correct or optimal way to imple- 
ment the VCA. For the case of local pseudopotentials, 
the implementation is straightforward the potential 
of the virtual system made from the (Ai- x Ti x )C alloy is 
generated simply by compositionally averaging the po- 
tentials of the parent AC and BC compounds, 

VVc-A(r) = (1 - x)V AC (r) + xV BC (r) . (1) 

In practice this is usually done in Fourier space by av- 
eraging Vac(G) and Vbc(G). In the case of semilo- 
cal (e.g, Hamann-Schluter-Chiang |ll|) pseudopotentials, 
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a similar averaging of the radial potentials Va,i(t) and 
VB,i{r) can be done separately in each angular mo- 
mentum channel I. However, with the fully non-local 
Kleinman-Bylander type separable pseudopotentials JTJ] 
that are most commonly used in the current generation 
of electronic-structure calculations, the implementation 
of the VCA is neither straightforward nor unique. For 
example, Ramer and Rappe fli"|| discuss four different 
ways of implementing the VCA for such pseudopoten- 
tials, each of them providing different physical results. 
Similarly, the best way of applying the VCA to the case 
of ultrasoft pseudopotentials |Q is less obvious still. 

The purpose of the present work is to report progress 
in addressing both of the above concerns. Taking them in 
reverse order, we first present a first-principles VCA ap- 
proach which is easily implemented for any type of pseu- 
dopotential. The method is demonstrated and tested in 
the context of calculations on Pb(Zr a; Tii_ a; )03 (PZT), 
an important perovskite solid solution. Our approach 
also has the advantage that it can easily be used to ob- 
tain estimates of the atomic forces that would arise if the 
real atoms were present, thus giving insight into the ex- 
pected displacements in the real alloy. Such information, 
which for example is highly relevant to many properties 
of ferroelectric systems, is not provided by the usual VCA 
techniques. 

Second, we use our new approach to evaluate the qual- 
ity of the VCA approximation for predicting the un- 
usual dielectric and piezoelectric properties in ferroelec- 
tric perovskite alloys. Perovskite compounds are known 
to exhibit anomalous dielectric properties. For exam- 
ple, they display anomalously large values of the Born 
effective charges, resulting from hybridization between 
the transition- metal d and oxygen 2p orbitals |]l6| , |l7| . 
Similarly, piezoelectric coefficients are large, compared to 
other classes of materials, both because of the large Born 
effective charges and because of the large microscopic re- 
action of the internal atomic coordinates to macroscopic 
strain @-|^). Using Pb(Zr 2; Tii_ 2; )03 with x=0.5 for 
our test system, Born effective charges and piezoelec- 
tric coefficients are calculated using our new VCA ap- 
proach together with the modern theory of polarization 
j2f| , p2f . We find that the VCA can be used with fair con- 
fidence to predict dielectric and piezoelectric properties 
of Pb(Zro.5Tio.5)03 alloys. As a matter of fact, our VCA 
technique yields large Born effective charges that are very 
nearly equal to the average between the effective charges 
of the parent compounds. Furthermore, comparison with 
calculations performed on ordered Pb(Zro. 5^0.5)03 su- 
percells demonstrates the ability of the VCA to mimic 
piezoelectric coefficients of perovskite solid solutions. 

The paper is organized as follows. In Sec. II, we im- 
plement our VCA approach in the context of density- 
functional theory, emphasizing the advantages of the new 
approach. Section III reports the predictions of this new 
VCA technique for the Born effective charges and piezo- 
electric coefficients of the Pb(Zr . 5^0.5)03 solid solution 
in its paraelectric and ferroelectric phases, respectively. 



We conclude in Sec. IV with a discussion of perspectives 
and future directions. The Appendix contains details 
about the implementation of our VCA technique within 
the Vanderbilt ultrasoft-pseudopotential scheme [p^| . 

II. THE VCA IMPLEMENTATION 

Within a pseudopotential approach to density-func- 
tional theory, the total energy of N v valence electrons 
can be written in terms of the one-particle wavefunctions 

(f>i as 

E tot [{4>i}, {Ri}] = C/({Ri}) + - v " + f«tl&> 

i 

4 + ^cN, (2) 

where 

y cxt (r, r') = £ V£(r - R,, r' - R,) , (3) 

R/ is the location of the site /, and Vp S are the pseu- 
dopotentials. Here, n(r) is the electron density, Exc is 
the exchange and correlation energy, and C/({Ri}) is the 
ion-ion interaction energy. A local pseudopotential takes 
the form y ext (r,r') = V cx t(v)5(r — r'), while a non-local 
pseudopotential is written as a sum of projectors. 

In either case, it is possible to derive a "VCA" operator 
equation by simply averaging the pseudopotentials of the 
alloyed elements on site /, 

V^(r,r') = (l-x)V p i(r,r')+xV p B s (r,r') , (4) 

where, e.g., A=Ti and B=Zr in Pb(Zr x Tii_ a: )03. For 
the sites occupied by the non-alloyed C elements [Pb or 
O in Pb^r^Tii-^Os], one simply takes 

V£{tJ)=V°(t,t>) . (5) 

Then V4 x t can be written 

Kxt(r,r')-^^<V7 s (r-R /Q ,r'-R /Q ) , (6) 

/ a 

where V^" is the pseudopotential for an atom of type a 
and is a "weight" which specifies the statistical com- 
position on site I. In cubic Pb(Zr x Tii_a;)03, for exam- 
ple, we set w = 1 for Pb and w — otherwise at the 
cube-corner site; w = I for oxygen and w = otherwise 
at the three face-center sites; and w — x and w = 1 — x 
respectively for Zr and Ti (and zero otherwise) at the 
cube-center site. In other words, we think of this crystal 
as composed of six atoms in the primitive cell: the usual 
one Pb and three oxygen atoms, and two "ghost" atoms 
(Zr and Ti) sharing the same lattice site (and having 
weights x and 1 — x respectively) . We then treat this unit 
cell containing six "atoms" in the usual first-principles 
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pseudopotential approach, solving the Kohn-Sham equa- 
tions in the presence of the potential given by Eq. (|6|) . 

This approach has the advantage of requiring only very 
slight modifications to the usual first-principles pseu- 
dopotential code. One simply inputs the weight w, along 
with the position and atom type, of each "atom" in the 
unit cell. These weights are then used in just a few places, 
e.g., in the construction of the total external potential 
. Perhaps the only subtlety is in the treatment of the 
Ewald energy [p3[ in the ion-ion interaction term U ap- 
pearing in Eq. (2). Here, we clearly have to prevent any 
Coulomb interaction between two "ghost atoms" on the 
same site, or else the Ewald energy would be infinite. 
In fact, the Ewald energy that we calculate is that of a 
crystal having valence charge 

a 

on site / (so that, in the case of isoelectronic substitution, 
this is just the usual Ewald energy). In practice, as long 
as the same gaussian width is used for all species when 
splitting the real-space and reciprocal-space Ewald con- 
tributions, this can be done very simply by replacing z" 
— > w a z" inside the program and deleting the infinite on- 
site interaction terms that would occur in the real-space 
Ewald sum. 

There are three definite advantages to this new VCA 
approach. First of all, it is extremely easy to imple- 
ment, as already indicated; only the minor modifications 
of Eqs. (^J) and (^) have to be implemented when start- 
ing from a conventional first-principles pseudopotential 
code. Second, in contrast to the approach of Ref. f p"l"| , 
there is no need to generate a pseudopotential for each 
virtual atom. Here, the pseudotentials arc created once 
and for all for each true atomic species; only the weights 
w/ change when dealing with a new composition of the 
solid solution. Third, the alloyed elements are still con- 
sidered as separate atomic species (with corresponding 
weights), rather than creating a single virtual atom as a 
whole. 

This last point is less trivial than it might seem. Be- 
cause the two "ghost atoms" on a site are considered as 
separate "atoms," one can consider responses to the dis- 
placement of just one of these "atoms" alone. In fact, 
the program automatically reports the forces on all the 
"atoms" in the unit cell, including those on the two ghost 
atoms separately. For structures of low symmetry, these 
forces need not be the same pj| . These forces can pro- 
vide some hints about the atomic distortions that would 
occur in the true disordered materials. The magnitudes 
of the force differences also act as a kind of internal diag- 
nostic for the appropriateness of using the VCA for the 
system of interest. A large magnitude would suggest that 
the VCA approximation is not expected to be very accu- 
rate, while a smaller value implies that the VCA can be 
used with fair confidence to mimic structural properties 
of the disordered alloy under consideration. 



It is also straightforward, using our approach, to use 
finite-difference methods to evaluate other kinds of re- 
sponse to the separate ghost-atom displacements. For 
example, the Born effective charges can be defined as 
the first-order polarization changes with respect to first- 
order sublattice displacements. As we will see in Sec. Ill, 
the contribution of each ghost atom to the Born effec- 
tive charge of the whole virtual atom can thus easily be 
calculated. To do the finite-difference calculation, we 
simply compute the change in polarization as the two 
ghost atoms are displaced to slightly different positions. 
Of course, to be meaningful, any real physical quantity 
(e.g., the derivative of some observable with respect to 
displacement) ultimately has to be evaluated at the con- 
figuration of identical ghost-atom positions. 

The discussion above has been limited to norm- 
conserving pseudopotentials, but the extension to the 
case of Vanderbilt ultrasoft pseudopotentials [[l5| is fairly 
straightforward. This extension is discussed in the Ap- 
pendix. In fact, all of our tests presented below have 
been carried out within the ultrasoft formulation. 



III. APPLICATION TO PEROVSKITE SOLID 
SOLUTIONS 

A. Born effective charges of the paraelectric 
Pb(Zr . 5 Tio. 5 )0 3 alloy 

Our first goal is to determine the dynamical effective 
charges of the Pb(Zro.soTio. 50)03 solid solution in its 
paraelectric phase, as predicted by the VCA. This alloy 
is usually denoted as PZT. For this purpose, we first per- 
form local-density approximation (LDA) [g5|] calculations 
within the Vanderbilt ultrasoft-pseudopotential scheme 
on the cubic perovskite structure, using our VCA tech- 
nique. As detailed in Ref. p6[ , a conjugate-gradient tech- 
nique is used to minimize the Kohn-Sham energy func- 
tional. The Pb 5d, Pb 6s, Pb 6p, Zr 4s, Zr Ap, Zr Ad, Zr 
5s, Ti 3s, Ti 3p, Ti 3d, Ti 4s, O 2s and O 2p electrons are 
treated as valence electrons. A weight w of 1 is assigned 
to Pb and oxygen atoms on their corresponding sites, 
while w = 0.5 for both Ti and Zr at the cube-center site 
[see Eqs. (6) and (7)]. Consequently, the VCA calcula- 
tion includes 44 electrons per cell. We use the Ceperley- 
Alder exchange and correlation p7{ as parameterized by 
Perdew and Zunger |jii8|| . A (6,6,6) Monkhorst-Pack mesh 
p9[ is used in order to provide converged results |2|] . The 
lattice parameter ao is fully optimized by minimizing the 
total energy, and is found to be very well described by 
Vegard's law. In other words, ao is very nearly equal to 
the compositional average between the lattice constants 
of pure PbTi0 3 and pure PbZr0 3 given in Ref. f||. 

To mimic the paraelectric phase, the atoms are kept in 
the ideal cubic positions. The dynamical effective charge 
Z33 a of each real and ghost atom is then calculated by 
using the formula 
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TABLE I. VCA dynamical effective charges Z 33 of 
Pb(Zro.5Tio.503), as compared with those of the parent com- 
pounds and their average. Last two rows show the comparison 
of the individual VCA Z 33 values for Zr and Ti with those of 
the parent compound [16]. 



Atom 


VCA 


PbZrOa 


PbTi0 3 


Ave 


Pb 


3.92 


3.92 


3.90 


3.91 


(B) 


6.47 


5.85 


7.06 


6.46 


Oi 


-2.54 


-2.48 


-2.56 


-2.52 


O3 


-5.28 


-4.81 


-5.83 


-5.32 


B: Zr 


9.62 


5.86 






B: Ti 


3.32 




7.06 





dP z = ^ Wa Z 33,a du z,a , (8) 

a 

where dP z is the change in polarization along the z- 
direction induced by the displacements du z<a of the a 
atoms along the z-direction, and w a refers to the weight 
assigned to the atoms. We allow the two ghost atoms 
to be at different atomic positions in order to compute 
the contribution of each of them on the Born effective 
charge of the whole virtual atom. We follow the pro- 
cedure introduced in Ref. plj which consists in directly 
calculating the spontaneous polarization as a Berry phase 
of the Bloch states. Technically, we use roughly 660 
Bloch states to assure a good convergence of the effec- 
tive charges. 

Table I shows the Born effective charges of the different 
atoms, as well as the compositional average between the 
Z33 of pure paraelectric PbTiC>3 and PbZrC>3 as given 
in Ref. In this Table, the averaged transition-metal 
atom interpolating between Zr and Ti is referred to as 
(£?), and the oxygen atoms are grouped into two kinds: 
those denoted O3, located between two (B) atoms along 
the z direction; and those denoted Oi, located between 
two (B) atoms in the perpendicular directions |L9|j30|| , 
One can first see that the VCA reproduces very well 
the effective charges of Pb, (B) and O atoms found in 
paraelectric PZT solid solutions Z| 3 is around 

4, —2.5, and —5.5 for Pb, Oi and O3 atoms, respec- 
tively, while the effective charge of the (B) atom is close 
to 6.5. In fact, Table I shows that the VCA approxi- 
mation essentially averages the effective charges of the 
parent compounds for Pb, Oi, O3 and (B). The VCA 
approach is thus able to mimic the weak and subtle in- 
teraction between the d orbitals of the transition-metal 
atom, treated as a single atom, and the O 2p orbitals. 
This interaction is responsible for the anomalous effec- 
tive charges of both O3 and (B) atoms [Tq] . On the 
other hand, Table I demonstrates that the Z* contribu- 
tions of the Zr or Ti ghost atoms in PbZro.5Tio.5O3 are 
not well approximated by the Z* 's of the corresponding 
Zr or Ti atom in the parent (PbZr03 and PbTiOs) com- 
pounds. This can be understood by realizing that the 
distance between Zr and O atoms in our VCA simula- 
tion of PbZro.5Tio.5O3 is smaller than the corresponding 



TABLE II. Structural parameters of Pb(Zr .5Tio. 5 03) 
within our VCA approach (denoted 'VCA') and for the super- 
cell ordered along the [100] direction used in Ref. [20] (denoted 
'Ordered 1'). Az are the ferroelectric atomic displacements, 
in c-lattice units. 





VCA 


Ordered 1 


Az(Pb) 


-0.0486 


-0.0480 


Az({B)) 


+0.0076 


+0.0064 


Az(Oi) 


+0.0790 


+0.0827 


Az(0 3 ) 


+0.0585 


+0.0555 



distance in PbZr03. The hybridization between the d 
orbitals of the Zr ghost atom and the O 2p orbitals then 
differs from the corresponding hybridization in PbZr03. 
This leads to an enhancement of the effective charge of 
the Zr ghost atom in PbZro.5Tio.5O3 with respect to the 
effective charge of Zr in PbZr03. Inversely, the effective 
charge of the Ti ghost atom in PbZro.5Tio.5O3 is much 
smaller than the effective charge of Ti in PbTiOs, since 
the distance between Ti and O atoms in the VCA calcu- 
lation of PbZro.5Tio.5O3 is larger than the corresponding 
distance in PbTi03. Interestingly, the underestimation 
for Ti cancels with the overestimation for Zr, yielding an 
effective charge of the whole (B) atom which is nearly 
equal to the average between the effective charges of the 
B atoms in the parents compounds. 

B. Piezoelectric coefficient of the ferroelectric 
Pb(Zr .sTio. 5 )0 3 alloy 

We now apply the VCA procedure to determine the 
piezoelectric coefficients e^ of the ferroelectric tetrago- 
nal P4mm ground state of the Pb(Zro.5oTio. 50)03 solid 
solution as predicted by the VCA. We use here the lat- 
tice constants minimiz ing the total energy of the or- 
dered supercell of Ref. pojj ; that is, the lattice parameter 
ao and the tetragonal axial ratio c/a are equal to 3.99 
A and 1.0345, respectively. The atomic displacements 
are fully optimized by minimizing the total energy and 
the Hellmann-Feynman forces, the latter being smaller 
than 0.05 eV/A at convergence. During these minimiza- 
tions, the ghost atoms can move but always share the 
same position. The Hellmann-Feynman force on the vir- 
tual (B) atom is simply the sum of the forces on the 
Zr and Ti "ghost" atoms, i.e., having a weight of 0.5 in 
Eqs. (6) and (7). In the ferroelectric VCA ground state 
of PZT, the force on the Ti atom is along the polariza- 
tion direction, i.e., along the z-axis. This force is exactly 
opposite to the force on the Zr atom. The magnitude of 
these forces is found to be 1.1 eV/A, which is less than 
three times larger than the force used to get convergent 
results in Ref. |3^]. This indicates that the VCA ap- 
proach can be used with some confidence to describe the 
properties of PZT alloys. As a matter of fact, Table II 
shows that the VCA can reproduce remarkably well the 
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TABLE III. Piezoelectric coefficients in C/m 2 of 
Pb(Zro. 5 Tio.50 3 ) within our VCA approach (denoted 'VCA'), 
for the supercell ordered along the [100] direction used in 
Ref. [20] (denoted 'Ordered 1'), for the supercell ordered along 
the [001] direction used in Ref. [18] (denoted 'Ordered 2'), 
and for the supercell ordered along the [111] direction used in 
Ref. [18] (denoted 'Ordered 3'). 





VCA 


Ordered 1 


Ordered 2 


Ordered 3 


e33 


4.4 


3.4 


4.8 


3.6 


e 3 3,c 


-0.8 


-0.8 


-0.7 


-0.7 


e33,i 


5.2 


4.2 


5.4 


4.3 



atomic displacements leading to the appearance of ferro- 
electricity in the ordered Pb(Zr .5oTi . 50)03 alloy 20|. 

Once the ferroelectric ground state is determined, the 
modern theory of polarization |^l],^2j is used to calculate 
the piezoelectric coefficient of PZT within our VCA pro- 
cedure. More precisely, the piezoelectric coefficients eij 
can be computed via fl34f 



1 

2^h 



^i? Q ,, — (ftG«-P) 



(9) 



where f2 is the cell volume and a — 1,2,3 runs over the 
three real-space lattice vectors R Q and reciprocal lattice 
vectors G Q , and r\j is the macroscopic strain. Eq. (9) has 
recently been derived in order to make the piezoelectric 
coefficients independent of the choice of branch of the 
Berry phase |34|. At the same time, Eq. (9) automati- 
cally eliminates of the so-called "improper" terms Jl9| as 
required to correctly predict the piezoelectric coefficients 
|p4| . Technically, Eq. (9) is evaluated by finite differ- 
ences between two strained configurations: first that of 
the ferroelectric ground state, and then for an additional 
1% strain relative to this ground state. In the second 
run, the relative atomic coordinates naturally have to be 
reoptimized in response to the applied strain. 

As done in Ref. |l^- po| , ^5| , the piezoelectric coefficients 
can be decomposed into "clamped-ion" and "internal- 
strain" contributions, 

e 33 = e 3 3, c + e 33ti 



(10) 



The "clamped-ion" or "homogeneous-strain" contribu- 
tion e33 iC is given by Eq. (9) for the case i = j = 3 and 
evaluated at vanishing internal strain (that is, without al- 
lowing the additional relaxation of the relative atomic co- 
ordinates that would be induced by the strain). e33 iC re- 
flects electronic effects, measuring the extent to which the 
Wannier centers fail to follow the homogeneous strain. 
The "internal-strain" part e 33i j measures just those con- 
tributions to the piezoelectric response coming from in- 
ternal distortions, i.e., reflecting the extent to which the 
ions fail to follow the homogeneous strain. In practice, 
633 and e33 jC are computed, and e3 3) , is then obtained 
from their difference. 



The results for 633, e33 jC and e 33) j, as predicted by the 
VCA, are shown in Table III and compared with various 
calculations on ordered supercells. One can first notice 
that the VCA is able to reproduce not only the magni- 
tudes but also the signs of the piezoelectric coefficients 
of ordered supercells. In particular, the clamped-ion con- 
tribution, which is negative and independent of the or- 
dering, is very well described by the VCA approxima- 
tion. Similarly, the VCA results for e 33i i and 633 lie be- 
tween those of the different ordered supercells, for which 
a larger spread exists. Overall, the results shown in Ta- 
bles II and III confirm the adequacy of the VCA to mimic 
ferroelectric properties of PZT. 



IV. CONCLUSIONS 

In summary, we have developed a new first-principles 
virtual crystal approach. This method (1) is easy to im- 
plement, (2) does not require the generation of pseudopo- 
tentials for each alloy composition, and (3) its outputs, 
via the computation of the Hellmann-Feynman forces on 
the "ghost" alloyed elements, provide a hint about the 
ability of the VCA to mimic properties of the disordered 
alloys under consideration. This technique has been 
applied, within the Vanderbilt ultrasoft-pseudopotential 
scheme Jl5| , to predict dielectric and piezoelectric prop- 
erties of the Pb(Zr . 5 Tio.5)03 solid solution in its para- 
electric and ferroelectric phase, respectively. Compar- 
ison with calculations performed on ordered supercells 
and with data on parent compounds demonstrates the 
adequacy of using the VCA for such properties, at least 
for isoelectronic perosvkite solid solutions. More work is 
needed to assess the ability of VCA to describe properties 
of heterovalent alloys pl] , p^ |, i.e., systems in which the al- 
loyed elements belong to different columns of the periodic 
table (e.g., Pb(Sco. 5 Nbo. 5 )03 or Pb(Mg 1/3 Nb 2/ 3)0 3 ). 

The present study also strongly suggests that first- 
principles derived effective-Hamiltonian methods, al- 
ready available for simple perovskite systems [^7|-^3[, 
can be used with confidence to predict finite-temperature 
properties of perovskite (isoelectronic) solid solutions, by 
modeling these alloys within the VCA approach. 



V. ACKNOWLEDGMENTS 

L.B. thanks the financial assistance provided by the 
Arkansas Science and Technology Authority (grant N99- 
B-21), and P. Thibado for the loan of a computer. 
D.V. acknowledges the financial support of Office of 
Naval Research grant N00014-97-1-0048. N.J. Ramer and 
A.M. Rappe are also thanked for communicating their re- 
sults |ll| to us prior to publication. Cray C90 computer 
time was provided at the NAVO MSRC under ONRDC 
# 1437. 



5 



VI. APPENDIX: IMPLEMENTATION OF THE 
VCA APPROACH WITHIN THE VANDERBILT 
ULTRASOFT-PSEUDOPOTENTIAL SCHEME 

In this appendix, we indicate how Eq. (6) can be re- 
alized when using the Vanderbilt's ultrasoft pseudo-po- 
tentials scheme flof . In this approach, the total energy 
of N v valence electrons described by the monoelectronic 
wavefunctions fa is given by 



E tot [{fa},{K I }} = ^2( 



Kxtl^i) 



dr dr 



, n(r)n(r') 



ExcM 



C/({Rx}) 



(11) 



where the external potential V cx t contains a local part 
V,Q° c n and a fully nonlocal part Vnl, 

V cxt = K° n + Vnl . (12) 
The local part contains local ionic contributions 

^°c n (r)-ECc U (|r-R7|) , (13) 
i 

while the fully nonlocal part is given by 

nm,I 

The functions as well as the coefficients 7 charac- 
terize the pseudopotentials, and thus differ for different 
atomic species. 

The electron density in Eq. (10) is given by 



nlr 



E 



rim, I 



Q I nm (r)(<Pi\0 I n )(P I m \<Pi) 



(15) 



where the augmentation functions Q^ m (r) are also pro- 
vided by the pseudopotentials and are strictly localized 
in the core regions. The ultrasoft pseudopotential is thus 
fully determined by the functions V™' 1 , Qn m and f3^, 
and by the scalar D^ m 7 . The algorithm used to gen- 
erate these quantities is described in Refs [^5[Q. The 
wavefunctions fa are eigensolutions of: 



H\fa) = eiS\fa) 
where S is an hcrmitian overlap operator given by 

s = i+ E iLKWll 



(16) 



nm.I 



with = J drQ* m (r), and where 

h = -v 2 + v cfi + £>LX></4 



Here V e g is the screened effective local potential 

VcS (r) = j^)= ^° n(r) + / dr 'y^v\ + ^ c(r) (19) 

with /i xc = 8Exc[n]/Sn(r), and 

Km = D^l + J dvVMQi^v) . (20) 

Eq. (6), which is the fundamental equation underlying 
our VCA approach, can thus be realized by simply re- 
placing three ionic quantities provided by the ultrasoft 
pseudopotentials, namely K* on ' a , Dnm, a , and Q" m , by 
their product with the corresponding atomic weight w a - 
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